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We combine large-scale atomistic modelling with continuum elastic theory to study the shapes 
of graphene sheets embedding nanoscale kirigami. Lattice segments are selectively removed from a 
flat graphene sheet and the structure is allowed to close and reconstruct by relaxing in the third 
dimension. The surface relaxation is limited by a nonzero bending modulus which produces a 
smoothly modulated landscape instead of the ridge-and-plateau motif found in macroscopic lattice 
kirigami. The resulting surface shapes and their interactions are well described by a new set of 
microscopic kirigami rules that resolve the competition between the bending and stretching energies. 
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Folding a two dimensional material lifts it into the 
third dimension enabling different physical functionali¬ 
ties. In a familiar example, folds can be introduced into a 
piece of paper to change its three dimensional shape with 
or without allowing for tears. A lattice model for the for¬ 
mer case (kirigami) has been studied recently [l! demon¬ 
strating rules for generating three dimensional shapes by 
the selective removal of segments from a parent honey¬ 
comb lattice and closing the tears by folding. Because 
the folding rules so defined are essentially geometrical it 
is possible that they could find applications in two di¬ 
mensional nanoscale materials and possibly even affect 
their electronic behavior [2]. 

In this Letter we examine this possibility by combining 
large scale atomistic modelling US] with analysis devel¬ 
oped from long wavelength elastic theory [SHU using a 
graphene sheet as a prototype. The models we adopt vi¬ 
olate two central tenets of macroscopic lattice kirigami: 
(a) the bending modulus is nonzero prohibiting the for¬ 
mation of sharply folded edges and (b) the medium is 
compressible allowing the system to store energy in shear 
and compressive strains. Thus, and perhaps not surpris¬ 
ingly, one finds that in a nanomaterial the regular faceted 
structures of macroscopic kirigami inevitably relax to 
softly rolling landscapes evocative of the English country¬ 
side. We find that these latter shapes, identified here in 
numerical calculations, can be accurately modelled and 
predicted using a new set of microscopic -kirigami rules 
appropriate in the weakly compressible limit. 

The shapes shown in Figure 1 compare two representa¬ 
tive kirigami-folded structures with their atomically re¬ 
laxed counterparts. The top panels contain defects in 
which atoms are removed from a strip and the gap is 
eliminated by rejoining lattice sites along a line that ter¬ 
minates on compensating edge dislocations containing 
nearest neighbor 5- and 7- membered rings. Fig. l(c- 
d) illustrate a deflection of this structure into the third 
dimensions via sharp folds that vertically displace the 
left and right hand regions in the same (panel (c)) or op¬ 


posite (d) directions. The defect energy density in this 
structure is confined to the edges of the folds so that the 
“up-up” ( uu ) and “up-down” (ud) patterns are degener¬ 
ate :lj. Starting from these structures we minimized the 
structural energy of a variety of atomistic models using 
interaction potentials for carbon developed by Los and 
Fasolino (LF) [3J0] which allow bonds to rupture and re¬ 
form and provide a useful description of the elastic prop¬ 
erties for carbon derived materials in diverse bonding en¬ 
vironments. The structures we develop should be con¬ 
trasted with patterned graphenes containing large open 
perforations designed to allow reversible large amplitude 
deformations under mechanical loading Hi- They are 
more akin to the fully bonded defect structures contain¬ 
ing height modulations found on scars that terminate on 
dislocation cores in some single layer graphenes produced 
by chemical vapor deposition nni- Two generic features 
of the fully relaxed structures are apparent in the lower 
panels of Figure l(e-f). For both defects we find a smooth 
variation in elevation that persists into the far field with 
soft pleats sourced by their near field defect structures. 

We quantify these observations by decomposing the 
height field h( r) on a disk of radius R into angular har¬ 
monics 

Mr) = X Mr)e im0 (1) 

m 

Figure 2(a) shows the radial dependence h m (r) for the 
allowed even m amplitudes in the shape in Figure 1(e). 
The relaxed structure is smooth, suppressing weight in 
its large m modes and confining its amplitude to the 
m = 0, ±2 deformations of the disk where h, 2 (r) (Fig. 
2(b)) is an increasing function of r out to the bound¬ 
ary. The bending energy has an areal energy density 
Ub = Kb(y 2 h ) 2 /2 and it is extremized by solutions of the 
biharmonic equation V 4 ft. = 0. We find that the radial 
dependences of our relaxed structures h m (r) are quite 
well described by linear combinations of these solutions 
projected into each angular harmonic subspace. For ex- 
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FIG. 1. Relaxation of a graphene kirigami in which atoms 
are removed from finite strip of flat graphene sheet (a) and 
the gap is rejoined, terminating on nearest neighbor pairs of 
5-7 disclinations (b). In (c) and (d) this structure is folded 
into three dimensions following the rules for macroscopic lat¬ 
tice kirigami where sharp folds separate plateaus that are dis¬ 
placed out of the plane in the same direction uu (c) or in op¬ 
posite directions ud (d). This microscopic structure relaxes to 
the shapes (e) and (f) generating a softly pleated landscape. 


FIG. 2. (a) The height field for the relaxed uu is decom¬ 

posed into its angular harmonics showing the radial depen¬ 
dence h m (r ) of its dominant contributions from m = 0, ±2. 
(b) The m = 2 radial dependence is well described by four 
solutions of a biharmonic equation projected into the m = 2 
subspace. The fit requires growing solutions with opposite 
signs which dominate the deflection in the far field. 


ample in the to = 2 subspace the representation 

h 2 (r) = h- 2 {r) = a 2 + + c 2 r 2 + d 2 r 4 (2) 

describes the shape as shown in Figure 2(c). Truncat¬ 
ing the expansion (1) to include only the m = 0 and 
to = ±2 solutions provides an excellent reconstruction 
of the exact shape as demonstrated in Figure 3(a). The 
ud structure (Fig. 1(f)) similarly relaxes to a smooth 
landscape well described by a superposition to = ±1, ±3 
angular harmonics. 

Note that the biharmonic equation admits two solu¬ 
tions that grow in the far field and generically these are 
both present in the relaxed structures but they always 
appear with opposite signs. Although it is tempting to 
attribute this to a boundary condition enforced at the 
edge of the disk, we find instead that this can be more 
easily understood as a global constraint on the shape. 
The growing solutions must compete in order to avoid 
a large strain energy penalty induced by their (locally) 
nonzero Gaussian curvatures. Note that a linear combi¬ 
nation of the growing solutions in Eqn. 2 make a contri¬ 
bution to the Gaussian curvature that is bilinear in the 
expansion coefficients for h 2 , explicitly we have for the 
determinant of the curvature tensor in the far field 

C 2 = —4 (02 + 6c 2 d 2 r 2 + &d\r A sin 2 (20)) (3) 

Following Nelson and Peliti [5] we recall that a coupling of 


the local Gaussian curvature to in-plane strain mediates 
nonlocal ultra-long range interactions between remote 
Gaussian curvatures, diverging in Fourier space oc q ~ 4 . 
Consequently, for a large system under open boundary 
conditions we can avoid a macroscopic energy that grows 
faster than the system size if its integrated Gaussian cur¬ 
vature vanishes. In the space of m-projected biharmonic 
solutions the residual Gaussian curvature cannot be made 
to vanish everywhere and with zero mean the residual 
curvature can be usefully described by its nonvanishing 
moments. For to = 2 and using Eqn. 3 we find that the 
disk-integrated curvature vanishes if the boundary ratio 
v = d 2 I{ 2 /c 2 = —0.423, in good agreement with the ratio 
(~ —0.47) obtained from our numerical calculations. We 
carried out similar analysis for different structures and in 
various angular momentum channels m in the expansion 
(1) and find that the boundary ratio is TO-dependent and 
consistent with our simulation data. 

The surfaces shown in Figure l(e,f) are therefore de¬ 
termined by three rules that resolve the competition be¬ 
tween its bending and stretching energy in the elasti¬ 
cally stiff (weakly compressible) limit: (1) the height 
field smooths by relaxing its amplitude to its low order 
symmetry-allowed angular harmonics, (2) the radial de¬ 
pendence in each ?n-channel superposes biharmonic solu¬ 
tions thereby producing a low bending energy, (3) these 
appear in “well-tempered” combinations that also avoid 
a large strain energy penalty by quenching the integrated 
Gaussian curvature. The defect energy is then deter- 
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mined by the core energy of the terminal dislocations, the 
bending energy in the extremal solution and the strain 
energy imposed by its residual Gaussian curvature. 

The argument given above fixes the amplitude ratio of 
the far field growing solutions but not their overall mag¬ 
nitudes which determines the degree of “warping” of the 
kirigami-ed disk. A scaling argument reveals that the lat¬ 
ter is determined by a boundary energy on the perimeter 
of the disc, presumably arising from the inequivalence 
of bulk (area) and surface (perimeter) interactions. For 
example, a structure with C 2 / 0 that results from a 
boundary interaction proportional to R and is opposed 
by a bulk interaction proportional to R? is described by 
an energy function 

U = aR 2 cl + /3 Rc 2 (4) 

where a > 0 and /3 are constants, giving c 2 = —/3/2aR. 
We can express the growing solutions of Eqn. 2 in a 
scaling form 



Thus for to = 2 by expressing all lengths (h, 7') in units of 
the disk radius R one obtains a universal warped shape 
determined by the value of /3. Note that this scaling rule 
is TO-dependent, i.e. different to’s all show scaling but are 
described by different scaling functions. The full shape is 
scalable to that the extent it can be described by a single 
dominant angular harmonic. In Figure 3(b) we test this 
hypothesis by plotting the scaled height h' = /R ver¬ 
sus the scaled radial coordinate r' = r/R demonstrating 
its near collapse to a single profile. We conclude that an 
unwarped kirigami profile with no growing solutions is 
nongeneric, and would require fine tuning the system to 
a special point at /? = 0. This is evidently not a property 
of the LF potentials for carbon [4] nor of any generic 
model for the interparticle interactions. Therefore the 
kirigami-ed disks generally feature a long distance shape 
modulation that cannot be confined to the defect. We 
interpret this as a microscopic analog to the step risers 
in macroscopic lattice kirigami that also propagate to 
the sample boundaries. It also suggests the possibility of 
tuning the shape of such a system by functionalizing the 
system boundaries as mechanism for controlling the edge 
potential parameter /3. 

These considerations can also be used to understand 
the energetics of microscopic kirigami. In macroscopic 
lattice kirigami the edges are sharp and the uu structure 
(Fig 1(c)) is degenerate in energy with the ud structure 
(Fig. 1(d)). Furthermore the energy of an uu config¬ 
uration is independent of the separation (d) of the dis¬ 
locations that define the vertices of their plateaus (Fig. 
1) since the sharp steps are nonoverlapping. These fea¬ 
tures do not apply to microscopic kirigami where the 
height profile is smooth and the dislocations can interact 
via overlap of their induced curvature fields. In Figure 



FIG. 3. (a) Reconstruction of the uu surface retaining only 

the m = 0, ±2 angular harmonics in the height field, (b) 
Numerical test of the scaling rule Eqn. 5 on four different 
disk radii. 


4(a) we compare the energies of the uu and ud configura¬ 
tions as a function of the vertex separation d. (To obtain 
these data the relaxation calculations were carried out on 
square rather than circular models so that the number of 
atoms is the same in each sampled structure.) The uu 
configuration is energetically preferred for any interver¬ 
tex spacing d, though for large d these energies converge 
to a common value which one can identify as twice the en¬ 
ergy of a single dislocation. At intermediate separations 
the energy degeneracy is in fact strongly broken, for ex¬ 
ample the energy difference for a separation of ~ 20A is 
« 0.5 eV. 

By analyzing these structures within continuum elas¬ 
tic theory we conclude that these energy differences arise 
from interactions that are mediated nearly entirely by the 
mean curvature of the extended overlapping height fields. 
The stretching energy, while present, is generally smaller 
than the energy stored in the mean curvature, and more 
importantly it is nearly s independent, indicating that 
its role is to simply renormalize the total self energy in 
these structures. The interactions between defects me¬ 
diated by the bending energy then lead qualitatively to 
the interaction profile shown in Figure 4. This behav¬ 
ior is captured even in a lowest order elastic theory. We 
first calculate the Lame coefficients A and and bend 
modulus Kb using our model potential giving the values 
presented in Table 1. In this expansion the energy can 
be partitioned into a pure bending contribution 

Ub = Y f d2v{ y 2h ) 2 ( 6 ) 
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and a strain term 

U a = 2 J + ^ u lk) (?) 

where Uij are the linearized strains (diUj +djui)/ 2. (We 
have investigated the role of the nonlinear strain terms 
that can appear Eqn. 7 and find that they do not qual¬ 
itatively change our conclusions.) Although the contri¬ 
bution from U s can be formally eliminated in favor of a 
(strongly) nonlocal interaction between Gaussian curva¬ 
tures 0, we choose instead to simply calculate the energy 
using the formula Eqn. 7. 


Elastic constant 

fitted value 

A 

3.23 eV/A 2 


10.67 eV/A 2 

B 

12.29 eV/A 2 

Kb 

93.49 eV 


TABLE I. Two dimensional Lame coefficients, bulk modu¬ 
lus and bending modulus obtained by fitting the structural 
energies for deformed graphene sheets using the interatomic 
potentials of Los and Fasolino [33 ■ 

In the continuum model one finds that energy degen¬ 
eracy of the uu and ud geometries is resolved and the 
uu configuration always favored. This can be under¬ 
stood if one regards the height fields of the two defects 
as additive. In the uu configuration the height deforma¬ 
tions appear with opposite signs and nearly cancel in the 
far field while in the ud configuration they interfere con¬ 
structively. The bending energy (though not the Gaus¬ 
sian curvature-induced stretching energy) is quadratic in 
derivatives of h and so the relative signs of the superposed 
height fields determines the sign of the their interaction. 
This behavior captures the essential results of the full 
atomistic calculations (Fig. 4). We also note that while 
the results obtained from the bare elastic theory correctly 
describe the ordering of the structural energies it fails to 
quantitively account for their magnitudes, as can be ex¬ 
pected since these structures are actually highly strained. 

Insights from the bending energetics of nanoscale 
kirigami may be useful for stabilizing structures in macro¬ 
scopic kirigami. The degeneracy of the uu and ud motifs 
is problematic for applications that would seek to stabi¬ 
lize a single target shape. This can be resolved by the 
introduction of macroscopic couplings that introduce an 
effective bending rigidity. Braces that suppress or pro¬ 
mote bending can be engineered to introduce nonlocal 
coupling between neighboring step risers and provide a 
route to encoding a unique surface structure. 

The analytic structure of our graphene-kirigami solu¬ 
tions also have important consequences for its Dirac elec¬ 
tronic structure near charge neutrality. In these struc¬ 
tures topological defects in their bond networks induce 



FIG. 4. Energies (expressed per unit area) for graphene 
kirigami as a function of vertex separation d. The uu and ud 
configurations are nondegenerate and the bend-induced ud 
potential is repulsive. These properties are described quali¬ 
tatively within a continuum elastic theory where the energy 
differences and their d dependence are controlled by the mean 
curvature in the relaxed structures. 


surface deformations with bend and (locally) nonzero 
Gaussian curvature. Separately, these structural features 
all couple to electronic motion in the tangent plane Hu¬ 
ng where the natural language for this coupling involves 
valley asymmetric bend- and strain-induced gauge fields 
EH. The gauge fields induced by pure bend are curl-free 
and have the innocuous effect of simply shifting the Dirac 
points in momentum space. By contrast Gaussian curva¬ 
ture is topologically nontrivial and links the system with 
a (valley dependent) local flux JOJi. The essential char¬ 
acteristic of the m-projected solutions presented above is 
that a competition between bending and stretching ener¬ 
gies generates a landscape where the Gaussian curvature 
is globally compensated (so that the total pseudo-flux is 
zero) but this can only be accomplished by sign changes 
on a network of nodal lines that carry the signature of 
the fully relaxed kirigami. The possibility of confining 
electronic modes along these lines and their role in defin¬ 
ing the low energy spectral and transport properties now 
presents an important problem for further study. 
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